In [ ]:

    

# Third-party
from astropy.io import ascii
import astropy.coordinates as coord
import astropy.units as u
import matplotlib as mpl
import matplotlib.pyplot as pl
import numpy as np
%matplotlib inline

import gary.coordinates as gc


    

In [ ]:

    

tbl = ascii.read("/Users/adrian/Downloads/xue08_bhb.txt")


    

In [ ]:

    

tbl.colnames


    

In [ ]:

    

c = coord.SkyCoord(l=tbl['GLON'], b=tbl['GLAT'], distance=tbl['d'], frame='galactic')
err_dist = u.Quantity(tbl['d']*0.1)


    

In [ ]:

    

vhel = u.Quantity(tbl['HRVel'])
err_vhel = u.Quantity(tbl['e_HRVel'])


    

In [ ]:

    

from ophiuchus import vcirc, vlsr, galactocentric_frame
from scipy.stats import norm
import emcee


    

In [ ]:

    

def ln_posterior(p, coordinate, obs_vlos, err_vlos, sigma_halo):
    pm_l,pm_b,vlos = p
    
    vgal = gc.vhel_to_gal(coordinate, 
                          pm=(pm_l*u.mas/u.yr,pm_b*u.mas/u.yr), 
                          rv=vlos*u.km/u.s,
                          vcirc=vcirc, vlsr=vlsr, galactocentric_frame=galactocentric_frame)
    vtot = np.sqrt(np.sum(vgal**2)).to(u.km/u.s).value
    return norm.logpdf(vtot, loc=0., scale=sigma_halo) + norm.logpdf(vlos, loc=obs_vlos, scale=err_vlos)


    

In [ ]:

    

# %%prun
# for i in range(100):
#     ln_posterior([1,1.,180.], c[0], vhel[0], err_vhel[0], sigma_halo)


    

In [ ]:

    

sigma_halo = 250. # km/s
ix = 0
sampler = emcee.EnsembleSampler(nwalkers=16, dim=3, lnpostfn=ln_posterior, 
                                args=(c[ix],vhel[ix].value,err_vhel[ix].value,sigma_halo))


    

In [ ]:

    

p0 = np.zeros((sampler.k, sampler.dim))
p0[:,0] = np.random.normal(0., 5., size=sampler.k) # mas/yr
p0[:,1] = np.random.normal(0., 5., size=sampler.k) # mas/yr
p0[:,2] = np.random.normal(vhel[ix], err_vhel[ix]/10., size=sampler.k) # km/s


    

In [ ]:

    

_ = sampler.run_mcmc(p0, 100)


    

In [ ]:

    

for i in range(3):
    pl.figure()
    pl.plot(sampler.chain[...,i].T, marker=None, drawstyle='steps', alpha=0.5)


    

In [ ]:

    

for i in range(3):
    pl.figure()
    pl.hist(np.hstack(sampler.chain[:,60:,i]))
    if i == 2:
        pl.axvline(vhel[0].value, c='#aaaaaa')
        pl.axvline(vhel[0].value+err_vhel[0].value, c='#aaaaaa', linestyle='dashed')
        pl.axvline(vhel[0].value-err_vhel[0].value, c='#aaaaaa', linestyle='dashed')


    

In [ ]:

    




    

In [ ]:

    

from scipy.integrate import tplquad


    

In [ ]:

    

def integrand(pm_l,pm_b,vlos,*args):
    p = (pm_l,pm_b,vlos)
    return np.exp(ln_posterior(p,*args))


    

In [ ]:

    

res = tplquad(integrand, 
              a=-500., b=500.,
              gfun=lambda *args: -100., hfun=lambda *args: 100.,
              qfun=lambda *args: -100., rfun=lambda *args: 100.,
              args=(c[0], vhel[0], err_vhel[0], sigma_halo))


    

In [ ]: